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Abstract — Point set registration is a key component in many computer vision tasks. The goal of point set registration is to assign 
correspondences between two sets of points and to recover the transformation that maps one point set to the other. Multiple factors, 
including an unknown non-rigid spatial transformation, large dimensionality of point set, noise and outliers, make the point set 
registration a challenging problem. We introduce a probabilistic method, called the Coherent Point Drift (CPD) algorithm, for both 
rigid and non-rigid point set registration. We consider the alignment of two point sets as a probability density estimation problem. We 
fit the GMM centroids (representing the first point set) to the data (the second point set) by maximizing the likelihood. We force the 
GMM centroids to move coherently as a group to preserve the topological structure of the point sets. In the rigid case, we impose the 
coherence constraint by re-parametrization of GMM centroid locations with rigid parameters and derive a closed form solution of the 
maximization step of the EM algorithm in arbitrary dimensions. In the non-rigid case, we impose the coherence constraint by regularizing 
the displacement field and using the variational calculus to derive the optimal transformation. We also introduce a fast algorithm that 
reduces the method computation complexity to linear. We test the CPD algorithm for both rigid and non-rigid transformations in the 
presence of noise, outliers and missing points, where CPD shows accurate results and outperforms current state-of-the-art methods. 

Index Terms — Registration, correspondence, matching, alignment, rigid, non-rigid, point sets, Coherent Point Drift (CPD), Gaussian 
mixture model (GMM), coherence, regularization, EM algorithm. 
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1 Introduction 

REGISTRATION of point sets is a key component in 
many computer vision tasks including stereo match- 
ing, content-based image retrieval, image registration 
and shape recognition. The goal of point set registration 
is to assign correspondences between two sets of points 
and/or to recover the transformation that maps one 
point set to the other. For example, in stereo matching, 
in order to recover depth and infer structure from a 
pair of stereo images, it is necessary to first define a 
set of points in each image and find the correspondence 
between them. An example of point set registration 
problem is shown in Fig. Q] The "points" in a point set 
are often features extracted from an image, such as the 
locations of corners, boundary points or salient regions. 
The points can represent both geometric and intensity 
characteristics of an image. 

Practical point set registration algorithms should have 
several desirable properties: (1) Ability to accurately 
model the transformation required to align the point sets 
with tractable computational complexity; (2) Ability to 
handle possibly high dimensionality of the point sets; 
(3) Robustness to degradations such as noise, outliers 
and missing points that occur due to imperfect image 
acquisition and feature extraction. 

The transformation usually falls into two categories: 
rigid or non-rigid. A rigid transformation allows only for 
translation, rotation and scaling. The simplest non-rigid 
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transformation is affine, which also allows anisotropic 
scaling and skews. Non-rigid transformation occurs in 
many real-world problems including deformable motion 
tracking, shape recognition and medical image regis- 
tration. The true underlying non-rigid transformation 
model is often unknown and challenging to model. 
Simplistic approximations of the true non-rigid trans- 
formation, including piece-wise affine and polynomial 
models, are often inadequate for correct alignment and 
can produce erroneous correspondences. Due to the usu- 
ally large number of transformation parameters, the non- 
rigid point sets registration methods tend to be sensitive 
to noise and outliers and are likely to converge into local 
minima. They also tend to have a high computational 
complexity. A practical non-rigid point set registration 
method should be able to accurately model the non-rigid 
transformation with tractable computational complexity. 

Multidimensional point sets are common in many real 
world problems. Most current rigid and non-rigid point 
sets registration algorithm are well suited for 2D and 
3D cases, but their generalization to higher dimensions 
are not always trivial. Furthermore, degradations such 
as noise, outliers and missing points significantly com- 
plicates the problem. Outliers are the points that are 
incorrectly extracted from the image; outliers have no 
correspondences in the other point set. Missing points 
are the features that are not found in the image due to 
occlusion or inaccurate feature extraction. A point set 
registration method should be robust to these degrada- 
tions. 

We present a robust probabilistic multidimensional 
point sets registration algorithm for both rigid and non- 
rigid transforms. We consider the alignment of two 
point sets as a probability density estimation problem, 
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Fig. 1 . The point set registration problem: Given two sets 
of points, assign the correspondences and the transfor- 
mation that maps one point set to the other. 

where one point set represents the Gaussian Mixture 
Model (GMM) centroids, and the other one represents 
the data points. We fit the GMM centroids to the data 
by maximizing the likelihood. At the optimum, the 
point sets become aligned and the correspondence is 
obtained using the posterior probabilities of the GMM 
components. Core to our method is to force GMM cen- 
troids to move coherently as a group, which preserves 
the topological structure of the point sets. We impose 
the coherence constraint by explicit re-parametrization 
of GMM centroid locations (for rigid and affine transfor- 
mations) or by regularization of the displacement field 
(for smooth non-rigid transformation). We also show 
how the computational complexity of the method can 
be reduced to linear, which makes it applicable for 
large data sets. The rest of the paper is organized as 
follows. In Section [2J we overview the current rigid and 
non-rigid point set registration methods and state our 
contributions. In Section [3j we formulate a probabilistic 
point set registration framework. In Sections |4] and |5j we 
describe our algorithms for rigid, affine and non-rigid 
registration cases, and relate them to existing works. In 
Section [6j we discuss the computational complexity of 
the method and introduce its fast implementation. In 
Section[7J we evaluate the performance of our algorithm. 
Section [8] concludes with some discussions. 

2 Previous Work 

Many algorithms exist for rigid and for non-rigid point 
set registration. They aim to recover the correspondence 
or the transformation required to align the point sets 
or both. Many algorithms involve a dual-step update, 
iteratively alternating between the correspondence and 
the transformation estimation. Here, we briefly overview 
the rigid and non-rigid point set registration methods 
and state our contributions. 

2.1 Rigid Point set Registration Methods 

Iterative Closest Point (ICP) algorithm, introduced by 
Besl and McKay [1] and Zhang [2], is the most popular 
method for rigid point set registration due to its simplic- 
ity and low computational complexity. ICP iteratively 
assigns correspondences based on a closest distance 
criterion and finds the least-squares rigid transformation 



relating the two point sets. The algorithm then redeter- 
mines the correspondences and continues until it reaches 
the local minimum. Many variants of ICP have been 
proposed that affect all phases of the algorithm from 
the selection and matching of points to the minimization 
strategy [3], [4]. ICP requires that the initial position of 
the two point sets be adequately close. 

To overcome the ICP limitations, many probabilistic 
methods were developed [5], [6]. These methods use 
soft-assignment of correspondences that establishes cor- 
respondences between all combinations of points accord- 
ing to some probability, which generalizes the binary 
assignment of correspondences in ICP. Among these 
methods are Robust Point Matching (RPM) algorithm 
introduced by Gold et al. [7], and its later variants [5], 
[8], [9]. In [10] it was shown that in RPM alternating 
soft-assignment of correspondences and transformation 
is equivalent to the Expectation Maximization (EM) al- 
gorithm for GMM, where one point sets is treated as 
GMM centroids with equal isotropic covariances and 
the other point set is treated as data points. In fact, 
several rigid point set methods, including Joshi and 
Lee [11], Wells [12], Cross and Hancock [13], Luo and 
Hancock [6], [14], McNeill and Vijayakumar [15] and 
Sofka et al. [16], explicitly formulate the point sets 
registration as a maximum likelihood (ML) estimation 
problem, to fit the GMM centroids to the data points. 
These methods re-parameterize GMM centroids by a 
set of rigid transformation parameters (translation and 
rotation). The EM algorithm used for optimization of the 
likelihood function consists of two steps: E-step to com- 
pute the probabilities, M-step to update the transforma- 
tion. Common to such probabilistic methods is the inclu- 
sion of an extra distribution term to account for outliers 
(e.g. large Gaussian [5] or uniform distribution [12]) and 
deterministic annealing on the Gaussian width to avoid 
poor local minima. These probabilistic methods perform 
better than conventional ICP, especially in presence of 
noise and outliers. 

Another class of rigid point sets registration methods 
are the spectral methods. Scott and Longuet-Higgins [17] 
introduced a non-iterative algorithm to associate points 
of two arbitrary patterns, exploiting some properties of 
Gaussian proximity matrix (Gram matrix) of point sets. 
The algorithm works well with translation, shearing and 
scaling deformations, but performs poorly with non- 
rigid transformations. Li and Hartley showed that corre- 
spondence and transformation are two factors of Gram 
matrices, and can be found iteratively using Newton- 
Schulz factorization [18]. This method performs well 
for moderate linear transformations. In spite of its el- 
egance, the large computational effort required for spec- 
tral methods prohibits its wide applicability There are 
a few other nonspectral methods worth mentioning. Ho 
et al. [19] proposed an elegant non-iterative algorithm 
for 2D affine registration by searching for the roots of 
the associated polynomials. Unfortunately this method 
does not generalize to higher dimensions. Belongie et 
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al. [20] introduced the "shape context" descriptor, which 
incorporates the neighborhood structure of the point set 
and thus helps to recover the correspondence between 
the point sets. 

Our approach to the rigid point sets registration is 
probabilistic and most closely related to the works of 
Rangarajan et al. [5], Wells [12] and Luo and Han- 
cock [14]. Despite extensive work in rigid probabilistic 
registration, none of the methods, to our best knowledge, 
provides a closed form solution to the maximization 
step (M-step) of the EM algorithm for a general mul- 
tidimensional case. The fact that the rotation matrix 
should be constrained to be orthogonal and to have a 
positive determinant further complicates its estimation. 
Rangarajan and collaborators [5] showed the solution for 
2D case only, where rotation is parametrized by a single 
angle. In higher dimensions the closed form solution 
with Euler angles parametrization is not feasible. Luo 
and Hancock [6], [14] find the rotation matrix through 
singular value decomposition (SVD). They ignored some 
terms of the objective function, which leads to only an 
approximate solution. We shall derive the exact closed 
form solution (M-step) for the rigid point set registration 
and discuss its difference from the related methods in 
Section |4] 

2.2 Non-rigid Point set Registration Methods 

Earlier works on non-rigid point set registration in- 
clude Hinton et al. [21], [22], who used the probabilistic 
GMM formulation. The GMM centroids are uniformly 
positioned along the contour (modeled using splines), 
which allows for non-rigid transformations. In practice, 
the method is applicable only to contour-like point 
sets. One of the most popular non-rigid point set reg- 
istration method is by Chui and Rangarajan [8], [9]. 
They proposed to use Thin Plate Spline (TPS) [23], [24] 
parametrization of the transformation, following RPM, 
which results into the TPS-RPM method. Similar to the 
rigid case, they use deterministic annealing and alternate 
updates for soft-assignment and TPS parameters estima- 
tion. They also showed that TPS-RPM is equivalent (with 
several modifications) to EM for GMM [10]. Tsin and 
Kanade [25] proposed a correlation-based approach to 
point set registration, which was later improved by Jian 
and Vemuri [26]. The method considers the registration 
as an alignment between two distributions, where each 
of the point sets represents the GMM centroids. One of 
the point sets is parametrized by rigid/ af fine parameters 
(in rigid /af fine case) or TPS (in non-rigid case). The 
transformation parameters are estimated to minimize the 
L2 norm between the distributions. These methods all 
use explicit TPS parametrization, which is equivalent 
to a regularization of second order derivatives of the 
transformation [23], [24]. The TPS parametrization does 
not exist when the dimension of points is higher than 
three, which limits the applicability of such methods. 

Huang et al. [27] proposed to implicitly embed the 
shape (or point sets in our case) into distance transform 



space, and align them similar to non-rigid image reg- 
istration algorithms. The authors use sum-of-squared- 
differences similarity measure between the embedded 
spaces and incremental free form deformation (FFD) to 
parameterize the transformation. The method performs 
well on relatively simple data sets. 

Finally, in our previous work we presented the Coher- 
ent Point Drift (CPD) algorithm [28] for non-rigid point 
sets registration. The algorithm regularizes the displace- 
ment (velocity) field between the point sets following the 
motion coherence theory (MCT) [29], [30]. Using varia- 
tional calculus, we obtained that the optimal displace- 
ment field has an elegant kernel form in multiple di- 
menstions. In this paper, we shall elaborate and analyze 
the CPD algorithm. We also extend the general non-rigid 
registration framework, and show that CPD and TPS- 
RPM are its special cases. Among other contributions, 
we estimate the width of GMM components without 
the time consuming deterministic annealing and show 
a fast CPD implementation to reduce the computational 
complexity to linear. We shall discuss and compare our 
method to the works of Chui and Rangarajan [9], and 
Jian and Vemuri [26] in Section [5] 

3 General Methodology 

We consider the alignment of two point sets as a proba- 
bility density estimation problem, where one point set 
represents the Gaussian mixture model (GMM) cen- 
troids, and the other one represents the data points. 
At the optimum, two point sets become aligned and 
the correspondence is obtained using the maximum of 
the GMM posterior probability for a given data point. 
Core to our method is to force GMM centroids to move 
coherently as a group to preserve the topological struc- 
ture of the point sets. Throughout the paper we use the 
following notations: 

• D - dimension of the point sets, 

• N, M - number of points in the point sets, 

• Xjv X £) = (xi, . . . , xjv) t - the first point set (the data 
points), 

• Y M xD = (yi, . . . ,y M ) T - the second point set (the 
GMM centroids), 

• T(Y, 6) - Transformation T applied to Y, where 6 
is a set of the transformation parameters, 

• I - identity matrix, 

• 1 - column vector of all ones, 

• d(a) - diagonal matrix formed from the vector a. 
We consider the points in Y as the GMM centroids, 
and the points in X as the data points generated by the 
GMM. The GMM probability density function is 

M+l 

p(x) = P{m)p(*\m) (1) 

m— 1 

where p(x|m) = ^ 27rrT l^ D / 2 exp ^ . We also added 
an additional uniform distribution p(x|M + 1) = to 
the mixture model to account for noise and outliers. We 
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use equal isotropic covariances a 2 and equal member- 
ship probabilities P(m) — jj for all GMM components 
(to = 1, ...,M). Denoting the weight of the uniform 
distribution as w, < w < 1, the mixture model takes 
the form 



1 M 1 

= w— + (1 - w) ^ JjP(*\ m ) 

m— 1 



(2) 



We re-parameterize the GMM centroid locations by a set 
of parameters 9 and estimate them by maximizing the 
likelihood, or equivalently by minimizing the negative 
log-likelihood function 



JV 



E(9,a 2 ) 



M+1 

^2 log ^2 p (™)p( x I to ) 



(3) 



where we make the i.i.d. data assumption. We define 
the correspondence probability between two points y m 
and x„ as the posterior probability of the GMM centroid 
given the data point: P(to|x„) = P(m)p(x n \m)/p(x n ). 

We use Expectation Maximization (EM) algorithm [31], 
[32] to find 6 and a 2 . The idea of EM is first to guess the 
values of parameters ("old" parameter values) and then 
use the Bayes' theorem to compute a posteriori proba- 
bility distributions P° ld (m\x. n ) of mixture components, 
which is the expectation or E-step of the algorithm. The 
"new" parameter values are then found by minimizing 
the expectation of the complete negative log-likelihood 
function [32] 

N M+1 

Q = ~Y,Y. P° ld (™\xn)log(P new (m)p n ™(x n \m)) 

n—1 m—1 

(4) 

with respect to the "new" parameters, which is called 
the maximization or M-step of the algorithm. The Q 
function, which we call the objective function, is also an 
upper bound of the negative log-likelihood function in 
((3). The EM algorithm proceeds by alternating between 
E- and M-steps until convergence. Ignoring the constants 
independent of 9 and a 2 , we rewrite (|Dl as 

N M 

Q(d,cr 2 ) = ^££ P ok V|x„) ||x„ -T(y m ,9)\\ 2 



N P D 



log a 2 (5) 



where JV P = £^=i £^=i P oM (to|x„) < JV (with JV = 
JVp only if w = 0) and P old denotes the posterior 
probabilities of GMM components calculated using the 
previous parameter values: 



1 II ^-T(y m ,»°") 



P oid (m|x„) 



exp 



E fc =i ex P 



(6) 



where c = (2na 2 ) D / 2 j^^-. Minimizing the function 
Q, we necessarily decrease the negative log-likelihood 
function E, unless it is already at a local minimum. To 



proceed, we specify the transformation T for rigid, affine 
and non-rigid point set registration cases separately. 

4 Rigid & Affine Point set Registration 

For rigid point set registration, we define the transforma- 
tion of the GMM centroid locations as T(y m ; R, t, s) = 
sRy m + t, where Rdxd is a rotation matrix, t£> x i is 
a translation vector and s is a scaling parameter. The 
objective function 0I takes the form: 

M,N 

Q(R,t, s , ( i 2 ) - — ^ P oM (to|x„) ||x„ - s Ry m - t|| 2 



Nt>D 



logo- 2 , s.t. R T R = I, det(R) = 1 (7) 

Note that the first term is similar to the one in the 
absolute orientation problem [33], [34], which is defined 

as minJ2n=i ll x « ~ ( s R-yn + 1)|| 2 in our notations. Equa- 
tion (O can be seen as a generalized weighted absolute 
orientation problem, because it includes weighted dif- 
ferences between all combinations of points. The exact 
minimization solution of the objective function (0 is 
complicated due to the constraints on R. To obtain the 
closed form solution we shall use Lemma [l] [35]: 

Lemma 1: Let Ru x d be an unknown rotation ma- 
trix and Afl X fl be a known real square matrix. Let 
US'S'V T be a Singular Value Decomposition (SVD) of 
A, where UU T = VV T = I and SS = d(s l ) with 
si > S2 >,...,> Sfl > 0. Then the optimal rotation 
matrix R that maximizes tr (A T R) is R = UCV T , where 
C = d(l,l,...,l,det(UV T )). 

To apply this lemma, we need to simplify the Q function 
to a form equivalent to tr(A T R). First, we eliminate 
translation t from Q. Taking partial derivative of Q with 
respect to t and equate it to zero, we obtain: 

t = -r^X T P T l - S R-^Y T P1 = /i x - sR/jy, 

iVp iVp 

where the matrix P has elements p mn = P old (m\x n ) in 
I0 and the mean vectors and /i y are defined as: 



E(X) 



— X T P T 1, u v = E(Y) = — Y T P1. 

JV fy N 



Substituting t back into the objective function and 
rewritting it in matrix form, we obtain 



Q = - T [tr(X T d(P T l)X) -2str(X T P T YR T ) 



s 2 tr(Y T d(Pl)Y) 



logcr 2 (8) 



where X = X — l/i x and Y = Y — 1/iy are the centered 
point set matrices. We use the fact that trace is invariant 
under cyclic matrix permutations and R is orthogonal. 
We can rewrite © as Q = ~a tr((X T P T Y) T R) + c 2 , 
where ci, c 2 are constants independent of R and c\ > 0. 
Thus minimization of Q with respect to R is equivalent 
to maximization of 

maxtr(A T R), A = X T P T Y, s.t. R T R = I, det(R) = 1. 
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Rigid point set registration algorithm: 

• Initialization: R = I, t = 0, s = 1, < w < 1 

2 _ 1 V^N II ||2 

G ~ DNM 2^n=l L^m=\ ll X ™ — ^ ™, \\ 

• EM optimization, repeat until convergence: 

• E-step: Compute P, 

_ 1 || x „-( s R,y m+ t)|| 2 

n - exp 2g 

ymn — T ii 7~Z ,,112 

^^ 1 C X p _ ^ l|X ' l_<SR ' yfc+t>11 + (2 7 T ( T 2 )D/ 2 _ML_M 

• M-step: Solve for R, s, t, a 2 : 
N P = l T Pl, Mx = iX T Fl )/ly = ^Y r Pl, 
X = X Y = Y - lfj%, 
A = X T P T Y, compute SVD of A = XJSSV T , 
R = UCV T , where C = d(l, .., 1, det(UV T )), 

_ tr(A r R) 
S ~~ tr(Y T d(Pl)Y)' 

t = /i X — SR/iy, 

a 2 = ^(tr(X r d(P T l)X) - str(A T R)). 

• The aligned point set is T (Y) = sYR T + lt T , 

• The probability of correspondence is given by P. 



Fig. 2. Rigid point set registration algorithm. 



Now we are ready to use Lemma [TJ and the optimal R 
is in the form 

R = UCV T , where U55V T = svd(± T P T Y) (9) 

and C = d(l, .., 1, det(UV T )). To solve for s and a 2 , we 
equate the corresponding partial derivative of (|8) to zero. 
Solving for R, s, t, a 2 is the M-step of the EM algorithm. 
We summarize the rigid point sets registration algorithm 
in Fig. 12 

The algorithm has one free paramater, w (0 < w < 1), 
which reflects our assumption on the amount of noise 
in the point sets. The solution for the rotation matrix is 
general D-dimensional. 

Affine point set registration: Affine registration case 
is simpler compared to the rigid case, because the 
optimization is unconstrained. Affine transformation is 
defined as T(y m ; R, t, s) = By m + t, where Rdxd is an 
affine transformation matrix, tu x i is translation vector. 
The objective function takes the form: 

M,N 

Q(B,t,„ 2 ) = — P old Mxn)||x„-(By m +t)|| 2 

m,n=l 

+ ^log. 2 (10) 

We can directly take the partial derivatives of Q, equate 
them to zero, and solve the resulting linear system of 
equations. The solution is straightforward and similar 
to the rigid case. We summarize the affine point set 
registration algorithm in Fig. [3] 

4.1 Related Rigid Point set Registration Methods 

Here, we discuss the probabilistic rigid point set regis- 
tration methods most closely related to ours. Rangarajan 
et al. [5] presented the RPM method for rigid point set 
registration. The method is shown for 2D case, where 



Affine point set registration algorithm: 

• Initialization: B = I,t = 0,0<w< 1 

2 _ 1 v^Af t-^ M || || 2 
a ~ DNM 2^n=l 2^m=l ll X ™ ~ 5^11 

• EM optimization, repeat until convergence: 

• E-step: Compute P, 

-^ll^-tBym+tJII 2 
n — c*P 2 " 

Fmn — T fi 77^ .,112 

££ll<=xp-5^ l|lc "- (Byfc+t)|1 +(27r CT 2 )«/ 2 T ^ 7 f 

• M-step: Solve for B,t,er 2 : 

• N P = l r Pl, Mx = J-X r P T l, % = ^Y^Pl, 

• X = X- Vj, Y = Y-l/iy 1 , 

• B = (X T P T Y)(Y T d(Pl)Y)- 1 , 

• t = /X X - B/iy, 

• a 2 = (tr(X T d(P T l)X) — tr(X T P T YB T )). 

• The aligned point set is T(Y) = YB T + lt T , 

• The probability of correspondence is given by P. 

Fig. 3. Affine point set registration algorithm. 

rotation matrix is parametrized by a single rotation 
angle, which allows to find an explicit update. Such Eu- 
ler 's angles approach is not feasible in multidimensional 
cases. RPM also uses deterministic annealing on a 2 , 
which requires to set the starting and stopping criteria 
as well as the annealing rate. The EM iterations has to be 
repeated at each annealing step, which can be slow. We 
argue that it is preferable to estimate a 2 instead of us- 
ing deterministic annealing. The deterministic annealing 
helps to overcome poor local minima, but for the rigid 
point set registration problem the rigid parametrization 
is a strong constraint that alleviates the advantages of 
the deterministic annealing. 

Luo and Hancock [14], [36] introduced the rigid point 
sets registration algorithm that is the most similar to 
ours. The authors optimize the objective function rather 
intuitively than rigorously, which leads to several as- 
sumptions and approximate minimization. They ignore a 
few terms of the objective function (see Eqs. 10,11 in [36]), 
where the last term does depend on transformation 
parameters, and must not be ignored. If such optimiza- 
tion converge, the M-step of the EM algorithm is only 
approximate. Among other differences, we want to men- 
tion that the authors use structural editing, a technique 
to remove some undesirable points, instead of using 
an additional uniform distribution to account for these 
points. Some other authors [15] also follow the rigid 
parameters estimation steps of Luo and Hancock [36]. 

5 Non-Rigid Point set Registration 

Non-rigid point set registration remains a challenging 
problem in computer vision. The transformation that 
aligns the point sets is assumed to be unknown and 
non-rigid, which is generally broad class of transforma- 
tions that can lead to an ill-posed problem. In order to 
deal with the problem we use Tikhonov regularization 
framework [37]-[39]. We define the transformation as the 
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initial position plus a displacement function v: 

T(Y,v)=Y + v(Y), 



(11) 



We regularize the norm of v to enforce the smoothness 
of the function [38]. Such approach is also supported by 
the Motion Coherence Theory (MCT) [29], [30], which 
states that points close to one another tend to move 
coherently, and thus, the displacement function between 
the point sets should be smooth. This is mathematically 
formulated as a regularization on the displacement (also 
called velocity) function. 

Additing a regularization term to the negative log- 
likelihood function we obtain 



(12) 



where E is the negative log-likelihood function 10, 4>(v) 
is a regularization term and A is a trade-off parame- 
ter. Such regularization is well formulated in Bayesian 
approach, where the regularization term comes from a 
prior on displacement field: p(v) = exp"*^"). 

We estimate the displacement function v using varia- 
tional calculus. We shall define the regularization term 
<fi(v) in different but equivalent forms and show that 
the optimal functional form of v is a linear combination 
of particular kernel functions. A particular choice of 
the regularization will lead to our non-rigid point set 
registration method. 

5.1 Regularization of the Displacement Function 

A norm of v in the Hilbert space H m is defined as: 



Hi 



E 

k=Q 



d k i 



dx k 



dx. 



(13) 



Alternatively, we can define the norm in the Reproduc- 
ing Kernel Hilbert Space (RKHS) [38], [40] as 



G(s) 



ds 



(14) 



where G is a unique kernel function associated with 
the RKHS, and G is its Fourier transform. Function v 
indicates the Fourier transform of the function v and s is 
a frequency domain variable. The Fourier domain norm 
definition has been used in the Regularization Theory 
(RT) [40] to regularize the smoothness of a function. 
Regularization theory defines smoothness as a measure 
of the "oscillatory" behavior of a function. Within the 
class of differentiable functions, one function is said to 
be smoother than another if it oscillates less; in other 
words, if it has less energy at high frequency. The high 
frequency content of a function can be measured by 
first high-pass filtering the function, and then measuring 
the resulting power. This can be represented by fl4)l , 
where G represents a symmetric positive definite low- 
pass filter, which approaches zero as ||s|| — > oo. For 
convenience, we shall write the regularization term as 



\Pv\\ 



(15) 



where an operator P "extracts" a part of the function for 
regularization, in our case, the high frequency content 
part [38], [39]. 

5.2 Variational Solution 

We find the functional form of v using calculus of varia- 
tion. Minimization of regularized negative log-likelihood 
function in (|T2l l boils down to minimization of the fol- 
lowing objective function (M-step): 



M,N 



Q(v, a 2 ) = JL P° ld (m\xn) K - (y m + v(y. 

N P D 



m,n—l 



logo- 



^\\PV\\ 2 



(16) 



A function v that minimizes (|T6l l must satisfy the Euler- 
Lagrange differential equation 



N M 

_££pold (m | Xn)(Xn _ (j 



v{y m )))S(z - y TO ) 
= PPv(z) (17) 



for all vectors z, where P is the adjoint operator to P. 
The solution to such partial differential equation can be 
written as the integral transformation of its left side with 
a Green's function G of the self-adjoint operator PP. 



M.N 



u ( z ) = £^° ld (™|x„)(x„-(y m +«(y m )))G(z,y m ) 



m.n— 1 



M 



= ^w m G(z,y m ) (18) 

m— 1 

where w m = ^ £n=i -P old (™|x n )(x„ - (y m +v(y m ))). 
Note that this solution is incomplete. In general, the 
solution also includes the term i/j(z) that lies in the null 
space of P [40], [41]. Thus, we achieve Lemma [2] 

Lemma 2: The optimal displacement function that 
minimizes < fl6] | is given by linear combination of the 
particular kernel functions centered at the points Y plus 
the term tp(z) in the null space of P: 



M 



v(z) = ^2 w m G(z,y m ) + ip(z) 



(19) 



m— 1 



where the kernel function is a Green's function of the 
self-adjoint operator PP. 

5.3 The Coherent Point Drift (CPD) Algorithm 

We choose the regularization term according to (14)1: 



ds 



(20) 



Ird G(s) 

where G is a Gaussian (note it is not related to the 
Gaussian form of the distribution chosen for the mix- 
ture model). There are several motivations for such a 
Gaussian choice: First, the Green's function (the kernel) 
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corresponding to the regularization term in d2"Dl is also 
a Gaussian (and remains a Gaussian for an arbitrary di- 
mensional case); the Gaussian kernel is positive definite 
and the null space term ip(z) = [40]. Second, by choos- 
ing an appropriately sized Gaussian function we have 
the flexibility to control the range of filtered frequencies 
and thus the amount of spatial smoothness. Third, the 
choice of the Gaussian makes our regularization term 
equivalent to the one in the Motion Coherence Theory 
(MCT) [30]: 



4>mct{v) 



E 

1=0 



l\2 l 



\D l v{x)\\ dx 



(21) 



where D is a derivative operator so that D 2l v = W 2l v and 
D 2l+1 v — V(V 2/ v), where V is the gradient operator and 
V 2 is the Laplacian operator. 

Lemma 3: The regularization term in l l20l with a 
Gaussian choice of low-pass filter G is equivalent to the 
the regularization term in 1 121) . Both terms represent the 
norm of the function v, after applying the operator P , 
and the corresponding Green's function is a Gaussian in 
both cases [38]. 

The equivalence of our regularization term with that 
of the Motion Coherence Theory implies that we are 
imposing motion coherence among the points and thus 
we call the non-rigid point set registration method the 
Coherent Point Drift (CPD) algorithm. 

We can obtain the coefficients w m by evaluating ((l"9t 
at y m points 

(G + Acr 2 d(Pl) _1 )W = d(Pl) _1 PX - Y (22) 

where Wmxd = (wi,...,wji/) t is a matrix of coeffi- 
cients, Gmxm is a kernel matrix with elements gij = 

i II y i~ y j II 2 

G(yi,yj) = e 2 I 13 II and d (•) is the inverse diag- 
onal matrix. The transformed position of y m are found 
according to ([TTJ as T = T(Y, W) = Y+GW. We obtain 
cr 2 by equating the corresponding derivative of Q to zero 

N M 



jT||x n -T(y m ,W) 



1 



N-pD 



N-pD 



(tr(X T d(P T l)X)-2 tr((PX) T T)+tr(T T d(Pl)T)) 



(23) 

We summarize the CPD non-rigid point set registration 
algorithm in Fig. [U 

Analysis: The CPD algorithm includes three free pa- 
rameters: w, A and (3. Parameter w (0 < w < 1) reflects 
our assumption on the amount of noise in the point 
sets. Parameters A and (3 both reflect the amount of 
smoothness regularization. A discussion on the differ- 
ence between A and (3 can be found in [29], [30]. Briefly 
speaking, parameter /3 defines the model of the smooth- 
ness regularizer (width of smoothing Gaussian filter in 
(|20] |). Parameter A represents the trade-off between the 
goodness of maximum likelihood fit and regularization. 



Non-rigid point set registration algorithm: 

^ M,N 

• Initialization: W = 0. a 2 = ||x„ 

DNM ^ " 

m,n— 1 

• Initialize w(0 < w < 1), (3 > 0, A > 0, 

- 1 II _ ||2 

• Construct G: g i3 = exp W* ny% y311 , 

• EM optimization, repeat until convergence: 

• E-step: Compute P, 

Pm.n — 



£f =1 exp 2„. 

• M-step: 



||*n-(yfc+G(*v)W)||- 



• Solve (G + Aa 2 ci(Pl)- 1 )W = d(Pl)- 1 PX - Y 
N P = 1 T P1, T = Y + GW, 

• a 2 = w i I7 (tr(X r d(P T l)X) - 2tr((PXf T)+ 

tr(T T d(Pl)T)), 

• The aligned point set is T = T(Y, W) = Y + GW, 

• The probability of correspondence is given by P. 



Fig. 4. The Coherent Point Drift algorithm for non-rigid 
point set registration. 



We note that solution of l|22|l gives the exact minimum 
of Q (16) , if cr 2 is assumed fixed. As far as we are 
estimating a 2 , l l2"2l and | |23) should be solved simulta- 
neously. The non-linear dependency of a 2 on W and 
vice-verse does not allow for simultaneous analytical 
solution. Iterative exact solution can be obtained by 
performing a few cyclic iterations on Ip2) and d23b within 
a single EM step. Practically, a single iteration, given by 
l|22) and |(23), decrease the Q function almost to the exact 
minimum. Such an iterative procedure, which decreases 
the Q function but not to exact minimum, is often called 
the generalized EM algorithm [31], [42]. 



5.4 Related Non-rigid Point set Registration Meth- 
ods 

The CPD algorithm follows our previous work [28] on 
non-rigid point sets registration. However, previously 
we have used deterministic annealing on a 2 , whereas 
here, we estimate the Gaussian width a 2 within ML 
framework. This allows us to significantly speed up 
the algorithm, alleviating the repeated EM-iterations 
for every single annealing step. We have not observed 
any decrease in accuracy of the method related to this 
change. In [28], we used a slightly different notation for 
the GMM centroid locations: we called Yo the initial 
centroids position (which we call Y here), and Y for 
the finall GMM centroid position (which we call T(Y) 
here). 

The most relevant non-rigid point sets registration 
algorithm to ours is TPS-PJ-TVI, more precisely its 
GMM formulation [10]. TPS-RPM uses Thin Plate Spline 
(TPS) [23], [24] parametrization of the transformation, 
which can be obtained by adding the regularization term 
that penalizes second order derivatives of the transfer- 



s 



mation. For instance, in 2D such regularization term is 




This term can be equivalently formulated in the Fourier 
space as: 

\\LTf= f ||s|| 4 |f(s)| 2 ds (25) 
Js. 2 

which is a special case of the Duchon splines [43]. The 
null space of such regularization includes affine trans- 
formations. Using the variational approach we can show 
that the optimal transformation T for such regularization 
is in the form T(Y) = YA + KC, where A is matrix of 
affine transformation cooefficients, C is a matrix of non- 
rigid cooefficients. For 2D case, matrix Kmx&t is kernel 
matrix with elements % = ||vj - yj || 2 log ||vj - Vj ■ ||. For 
3D case, matrix K has elements fcy = ||y^ — [[. For 
4D or higher dimensions the TPS kernel solution does 
not exist [44]. Finally, to link such regularization to 
our non-rigid registration framework, we note that the 
regularization of the displacement field v, instead of the 
transformation itself, is exactly the same, because, (|24l 
is invariant under affine transformations, in other words 
\\LT(z)\\ 2 = ||L(z + v(z))|| 2 = \\Lv(z)\\ 2 . This means that 
both CPD and TPS-RPM regularizes the displacement 
function, but using different regularization terms. 

The advantage of CPD regularization (as given by ((20) 
or ll2il ) comparing to TPS ( (pit or l|25|l ), is that it easily 
generalizes to N dimensions. Also we can control the 
locality of spatial smoothness by changing the Gaussian 
filter width (3, whereas TPS does not have such flexibility. 
This, however, introduces one extra parameter to the 
method, but TPS-RPM has to uses one extra parameter 
to regularize affine matrix after all. Among other differ- 
ences, TPS-PvPM approximates the M-step solution of the 
EM algorithm [10] for simplicity and use deterministic 
annealing on a 2 . 

Finally, Jian and Vemuri [26] consider the registration 
as an alignment between the distributions of two point 
sets, where a separate GMMs are used to model the 
distribution for the point sets. One of the point sets 
is parametrized by TPS. The transformation parameters 
are estimated to minimize the L2 norm between the 
distributions. In our case, the CPD method maximizes 
the likelihood function, which is equivalent to KL diver- 
gence minimization between two mixture distributions: 
GMM and mixture of delta functions. KL divergence is 
more appropriate similarity measure for the densities 
than L2 norm, because it weights the error according 
to its probability. 

6 Fast Implementation 

Here we show that CPD computational complexity can 
be reduced to linear up to a multiplicative constant. 
We use the fast Gauss transform (FGT) [45] to compute 
the matrix-vector products PI, P T 1, PX, which are the 
bottlenecks for both rigid and non-rigid cases. We use 



Compute P T 1,P1 and PX: 

• Compute K T 1 (using FGT), 

• a = l./(K T l + cl), 

• P T 1 = 1 - ca, 

• PI = Ka (using FGT), 

• PX = K(a. * X) (using FGT), 

Fig. 5. Matrix-vector products computation through FGT. 

low-rank matrix approximation to speed-up the solution 
of the linear system of equations d22l l for the non-rigid 
case. 

The fast Gauss transform: Greengard and Strain [45] 
introduced the fast Gauss transform (FGT) for fast com- 
putation of the sum of exponentials: 

N 

f(y m ) = ^exp-^ l|x "- ym|12 , Vy m , m=l,...,M. 

n=l 

The naive approach takes O(MN) operations, while FGT 
takes only 0(AI + N). The basic idea of FGT is to expand 
the Gaussians in terms of truncated Hermit expansion, 
and approximate 10 up to the predefined accuracy. 
Rewriting I0 in matrix form, we obtain f = Kz, where z 
is some vector and K.mxn is a Gaussian affinity matrix 
with elements: k mn — exp~2^" x "~ r< ' y " 1 ' ) " , which we 
have already used in our notations. We simplify the 
matrix-vector products PI, P T 1 and PX, to the form 
of Kz and apply FGT. Matrix P © can be partitioned 
into 

P = Kd(a), a= l./(K T l + cl) (26) 

where d(a) is diagonal matrix with a vector a along the 
diagonal. Here, we use Matlab element-wise division (./) 
and element-wise (.*) multiplication notations. We show 
the algorithm to compute the bottleneck matrix-vector 
products PI, P T 1 and PX using FGT in Fig. [5] We note 
that for dimensions higher than three, we can use the im- 
proved fast Gauss transform (IFGT) method [46], which 
is a faster alternative to FGT for higher dimensions. 

During the finall EM iterations, the width of the Gaus- 
sians a 2 becomes small. The Hermitian expansion thus 
requires many terms to approximate highly multimodal 
Gaussian distribution for a given precision. At the final 
iterations, the Gaussian becomes very narrow, and we 
can switch to the truncated Gaussian approximation (set 
zeros outside a predefined box). 

Low-rank matrix approximation: In the non-rigid case, 
we need to solve the linear system (22), which is 0(M 3 ) 
using direct matrix inversion. We note that the left hand 
side matrix of | |22|| is symmetric and positive definite. 

We use low-rank matrix approximation of G, where 
G is a Gaussian affinity matrix with elements = 
exp~2^" yi ~ yj " . We approximate the matrix G as 

G = QAQ T (27) 

where ^k-kk is a diagonal matrix with K largest eigen- 
values and the matrix Qmxk is formed from the cor- 
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C) 



Initialization Iteration 10 Iteration 30 Iteration 40 Result (iteration 50) 

Fig. 6. Fish data set, rigid registration examples. We align Y (blue circles) onto X (red stars). The columns show 
the iterative alignment progress, a) Registration of the point sets with missing non-overlapping parts (w = 0.5); b) 
Registration of the point sets corrupted by random outliers (w = 0.5); c) A challenging rigid registration example, 
where both point sets are corrupted by outliers and biased to different sides of the point sets. We have also deleted 
some parts from both point sets. We set w = 0.8 and fix scaling s = 1. CPD registration is robust and accurate in all 
experiments. 



responding eigenvectors. G is the closest AT-rank ma- 
trix approximation to G both in L2 and Frobenius 
norms [47]. To solve the linear system in j22|l we use 
the Woodbury identity and invert the first term as 

(QAQ T + Xa 2 d(Pl)" 1 )- 1 = t\ d(Pl) 

- J^y (i(Pl)Q(A- 1 + ^Q T d(Pl)Q)- 1 Q T d(Pl) 

(28) 

The inside matrix inversion is of 0(AT 3 ), where K <c 
M. For instance choosing K = A/ 1 / 3 largest eigenvalues, 
we reduce the computational complexity to linear. We 
can pre-compute K largest eigenvalues and eigenvectors 
of G using deflation techniques [48]. It requires several 
iterations with the matrix-vector product Gz, which can 
be implemented explicitly or through FGT. 

The low-rank matrix approximation intuitively con- 
straints the space of the non-rigid transformations, and 
can be even desirable to further constrain the non-rigid 
transformation. If the number of points is large and well 
clustered, then an extremely small percent of eigenvalues 
will be sufficient for an accurate approximation. 

7 Results 

We implemented the algorithm in Matlab, and tested it 
on a Pentium4 CPU 3GHz with 4GB RAM. We imple- 
mented the matrix-vector products in C as a Matlab mex 



functions to avoid the storage of P. The code is avail- 
able at www.csee.ogi.edu/~myron/matlab/cpd. 
We shall refer to our method as Coherent Point Drift 
(CPD) both for rigid and non-rigid point sets registration 
methods presented in this paper. We have also imple- 
mented the matrix-vector products through FGT using 
the Matlab FGT implementation by Sebastien Paris [49]. 

We consider rigid and non-rigid experiments sepa- 
rately below. We shall always pre-align both point sets 
to zero mean and unit variance before the registration. 

7.1 Rigid Registration Results 

We show the CPD rigid registration on several examples, 
test the fast CPD implementation and evaluate its per- 
formance in comparison with LM-ICP [3], which is one 
of the most popular robust rigid point set registration 
methods. 

Rigid fish point set registration: Fig. [6] shows several 
rigid regsitration tests on 2D fish point sets. In Fig. [6^ 
we deleted non-overlapping parts in both point sets 
and set w = 0.5, where w is a weight of the uniform 
ditribution that accounts for noise and outliers. In Fig. [6j) 
we corrupted the point sets by outliers. We generate 
outliers randomly from a normal zero-mean distribution. 
CPD demonstrates robust and accurate performance in 
all examples. Fig. [6fc demonstrates a challenging exam- 
ple, where both point sets have missing points and are 
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Initialization Iteration 10 Iteration 20 Iteration 30 Result (iteration 50) 



Fig. 7. 3D bunny point set rigid registration examples. We align Y (blue circles) onto X (red dots). The columns show 
the iterative alignment progress. We initialized one of the point sets with 50 degree rotation and scaling equal 2. a) 
Registration of the point sets with missing points (w = 0.5); b) A challenging example of CPD rigid registration with 
missing points, outliers and noise. CPD shows robust and accurate registration result in all experiments. 




Performance 0.04 Noise STD 0.12 Noise STD 0.2 Noise STD 



Fig. 8. A comparison of CPD and LM-ICP rigid registration performances with respect to noise in the X (first row) and 
the Y point sets (second row). We align Y (blue circles) onto X (red dots). The columns 2,3 and 4 show the examples 
of initial point sets for different random noise stds added to the point set positions. The first column shows the error in 
estimating the rotation matrix for CPD (blue) and TPS-RPM (red). CPD outperforms LM-ICP in all cases. 



corrupted by outliers. The most challenging here is that 
we biased the outliers to the different sides of fish point 
sets. We were able to register such point sets only by 
fixing the scaling to be constant (estimating rotation and 
translation only). CPD demonstrates accurate and robust 
registration performance. 

Rigid bunny point set registration: We test 3D rigid 
point sets registration on the Stanford "bunny" data 
set [50]. We use a subsampled bunny version of 1889 x 3 
points. In Fig. [7^, we have deleted the front and back 
parts of the bunny point sets. In Fig. 03, we have added 
random outliers to different sides of the point sets. We 
set w = 0.7. CPD registration is accurate and robust in 
all examples. 

We compare the CPD rigid algorithm to the LM-ICP 
method [3], a robust version of ICP Fig. [8] shows the 
performance of CPD and LM-ICP with respect to noise 
in the point sets. We align the Y point set (blue circles) 
onto the X point set (red dots). We set w = 0.5. The 
known initial rotation discrepancy between the point 



sets is 50 degrees. The first and second rows shows the 
alignment performance when a random noise is added 
to the X and Y point set positions respectively. We use 
a norm of the difference between the true and estimated 
rotation matrix as an error measure. A few initial point 
sets examples with different noise std are shown in the 
columns 2, 3 and 4 of Fig. [8] For each level of the noise 
stds we made 25 independent runs. The first column 
plots the error values (mean and standard deviation) 
in the estimated rotation matrix as a function of noise 
levels. The CPD rigid algorithm outperforms the robust 
LM-ICP method, especially when the noise is present in 
the X point set. 

Fig- shows the performance of CPD and LM-ICP 
with respect to the outliers in the point sets. We add 
different number of outliers (irrelevant random points) 
to the point sets. An examples of such initial point sets 
are shown in columns 2, 3 and 4 of Fig. [9] for 600, 1800 
and 3000 outlier points added respectively. The first and 
second row show the cases of outliers present in the 
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Performance Outliers 600 Outliers 1800 Outliers 3000 



Fig. 9. A comparison of CPD and LM-ICP rigid registration performances with respect to outliers in the X (first row) 
and the Y (second row) point sets. We align Y (blue circles) onto X (red dots). The columns 2,3 and 4 show the 
examples of initial point sets with different number of outliers added. The first column show the error in estimating the 
rotation matrix. CPD outperforms LM-ICP. 




c) 



Initialization Iteration 10 Iteration 20 Iteration 40 Result (iteration 50) 

Fig. 10. Non-rigid CPD registration of 2D fish point sets, a) Noiseless fish point sets registration (91 x 2 points, w = 0); 
b) Registration of 2D fish point set with missing points (w = 0.5); c) Registration of 2D fish point set in presence of 
outliers (w = 0.5). CPD registration is robust and accurate in all experiments. 



X and Y point sets respectively. CPD performs well 
in all experiments, whereas LM-ICP performance is less 
accurate. 

Fast rigid CPD implementation: We also test the CPD 
performance with FGT implementation of the matrix- 
vector products. We use four Stanford bunny sets of 
sizes: 453 x 3, 1889 x 3, 8171 x 3 and 35947 x 3. For each 
of the cases we add a small amount of noise and outliers 
to both point sets, initialized them with 50 degrees 
rotation and set w — 0.3. For the FGT parameters, 
we used "ratio of far field"=8, "number of centers"=80, 



"order of truncation"=5. Table. [TJ shows the registration 
time with and without FGT. The FGT implementation 
is significantly faster. We note that there are several 
downsides of using the FGT: a) FGT requires its own 
parameter initialization; b) CPD (with FGT) aligns the 
point sets to 0.1 degree error rotation and then starts 
being unstable. This is because a 2 becomes small and 
the FGT approximation error becomes significant. At 
this point one can either stop (the alignment already is 
reasonably accurate) or proceed with ICP or truncated 
Gaussian CPD. 
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Deformation level 




Fig 



(a) (b) (c) 

1 1 . A comparison of CPD and TPS-RPM on the 2D fish point sets with respect to a) Deformation level; b) Noise 



level; c) Outliers. CPD shows more accurate registration performance compared to TPS-RPM, especially in presence 
of outliers and complex non-rigid deformations. 



N, M 


Naive 


FGT 


453 x 3 


0.6s 


0.7s 


1889 x 3 


lis 


3s 


8171 x 3 


4m 


10s 


35947 x 3 


3.5hr 


51s 



TABLE 1 

The rigid CPD registration time for naive (no FGT) and 
FGT implementations. The FGT-based implementation is 
significantly faster. 



N,M x D 


Naive 


FGT 


Low-rank 


FGT & Low-rank 


453 x 3 


2s 


2.3s 


1.7s 


1.8s 


1889 x 3 


lm22s 


1ml 6s 


19s 


lis 


8171 x 3 


3hr 


2hr26m 


10m20s 


lm37s 


35947 x 3 






40m 


10m 



TABLE 2 

Registration time required for non-rigid registration of 3D 
bunny point sets. The time is shown when using only 
FGT of vector-matrix products, only low-rank matrix 
approximation of Gaussian kernel matrix or both. 



7.2 Non-rigid Registration Results 

We show CPD non-rigid registration on several exam- 
ples, test the fast CPD implementation and evaluate CPD 
performance in comparison to TPS-RPM [9], which is one 
of the best performing non-rigid point set registration 
methods. We set A = 2, (3 = 2. 

Non-rigid fish point set registration: Fig. [TOa shows 
non-rigid CPD registration of two fish point sets with 
clean data. Fig. [lOb is with missing points (w = 0.5). 
Fig. [TOc is with both point sets are corrupted by outliers 
(w — 0.5). The non-rigid CPD registration results are 
accurate in all experiments. 

We test CPD against TPS-RPM [9] on synthetically 
generated 2D fish non-rigid examples with respect to 
a) level of non-rigid deformation, b) amount of noise 
in the point sets locations c) number of outliers. We 
set w — 0.3 in all experiments. Since we know the 
true correspondences, we use the mean squared distance 
between the corresponding points after the registration 
as an error measure. For each set of parameters we 
have conducted 25 runs. Fig. [TTh shows the methods 
performances with respect to the level of initial non- 
rigid deformation between the point sets. To generate 
the non-rigid transformation, we parameterize the point 
sets domain by a mesh of control points, perturb the 
points and use splines to interpolate the deformation. 
The higher level of mesh point perturbations produce 
the higher deformation. CPD shows accurate registration 
performance and outperforms the TPS-RPM. Fig. fTTb 
shows the methods performances with respect to the 
amount of noise. We add a zero-mean white noise with 
increasing levels of stds to the point sets. Both CPD and 



TPS-RPM show accurate performances. We note that, 
due to deterministic annealing used by TPS-RPM, its 
convergence takes significantly more iterations and time. 
Fig. [TTh shows the methods performances with respect 
to the number of outliers. We add random outliers to the 
point sets and plot the registration error with respect to 
the ration of number of outliers to the number of data 
points. CPD shows robust registration performance and 
outperforms the TPS-RPM. 

Non-rigid 3D face registration: We show the CPD per- 
formance on 3D face point sets. Fig. [12h shows two 3D 
face point sets related through non-rigid deformation. 
Fig. [12b shows two 3D face point sets point sets with 
added outliers and non-rigid deformation. Non-rigid 
CPD registration is accurate in all experiments. 
Non-rigid 3D LV point set registration: Finally, we 
demonstrate the CPD performance on non-rigid a 3D left 
ventricle (LV) contours segmented from 3D ultrasound 
images, using active contour based segmentation [51]. 
Fig. [13] shows (a) two LV point sets at different time 
instances, (b) the registration result, (d) the displacement 
field required for CPD alignment. That the registration 
result is accurate. 

Fast non-rigid CPD implementation: We test the com- 
putational time of the fast CPD non-rigid implementa- 
tion on several subsampled 3D Stanford bunny point 
sets. We use FGT of the matrix-vector products, the 
low-rank matrix approximations of the kernel matrix, or 
both. We applied a moderate non-rigid deformation to 
the bunny point sets. The registration time of the non- 
rigid CPD is shown in Table [2] We were unable to run 
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a) 




Initialization Result 

Fig. 12. Non-rigid registration of 3D face point sets, a) 
Registration of clean point sets b) Registration of point 
set with outliers. CPD shows accurate alignment. 




(a) Initialization (b) Result (c) Displacement 

Fig. 13. Non-rigid registration of 3D left ventricle (LV) 
point sets, (a) two LV point sets at different time instances, 
(b) the registration result, (c) displacement field between 
the corresponding points. 




Fig. 14. The log-plot of the eigenspectrum of the kernel 
matrix G for the bunny point sets of size 1889 x 3. 



the test without the low-rank matrix approximation for 
the largest bunny set (35947 x 3), because of the large 
RAM requirements to construct the kernel matrix G. We 
used only 100 leading eigenvalues and eigenvectors in 
all cases. Table|2]shows that the main computational bot- 
tleneck is in solving the linear system of equations l(22)l . 
because the low-rank matrix approximation alone can re- 
duce the computational time significantly Both FGT and 
low-rank approximations provide further speed-up with 
only moderate loss of accuracy. We note that almost 60% 
of the time required to complete the CPD registration 
using the low-rank matrix approximation were required 
to pre-compute the eigenvalues and eigenvectors of the 



kernel matrix G. 

We also show the eigenvalues for a particular example 
of the bunny point set of size 1889 x 3 in Fig. [14] 
Eigenvalues drops quickly below 10~ 3 only after 10 
largest eigenvalues, and drops below 10~ 5 after about 
100 eigenvalues. The approximation error of using a low 
rank approximate matrix (constructed from 100 leading 
eigenvectors and eigenvalues), is only 10~ 8 in terms of 
its Frobenius norm. 

8 Discussion and Conclusion 

We introduce a probabilistic method for rigid and non- 
rigid point set registration, called the Coherent Point 
Drift algorithm. We consider the alignment of two point 
sets as a probability density estimation, where one point 
set represents the Gaussian Mixture Model centroids, 
and the other represents the data points. We iteratively 
fit the GMM centroids by maximizing the likelihood and 
find the posterior probabilities of centroids, which pro- 
vide the correspondence probability. Core to our method 
is to force the GMM centroids to move coherently as a 
group, which preserves the topological structure of the 
point sets. 

Our contribution includes the following aspects. For 
the non-rigid point set registration, we formulate the 
motion coherence constraint and derive a solution of 
the regularized ML estimation through the variational 
approach, which leads to an elegant kernel form. CPD 
simultaneously finds both the transformation and the 
correspondence between two point sets without making 
any prior assumption of the non-rigid transformation 
model except that of motion coherence. For the rigid 
case, we derived a closed form multidimensional solu- 
tion (of the M-step of the EM algorithm), which has not 
been derived exactly before. Finally, we introduced the 
fast CPD implementation using fast Gauss transform and 
low-rank matrix approximation to reduce the computa- 
tional complexity of the method to as low as linear. On 
top of the computational advantage, the low-rank kernel 
approximation provides more stable solutions in cases 
where the matrix G is poorly conditioned. To our best 
knowledge, CPD is the only method capable of non-rigid 
registration of large data sets. Both rigid and non-rigid 
CPD registration methods are multidimensional and can 
be applied to arbitrary dimensional data sets. 

We estimate the GMM width, a 2 , within the ML 
formulation. We have not observed any decrease in 
performance compared to the deterministic annealing 
approach. Estimation a 2 allows to reduce the number 
of free parameters and, most importantly, to significantly 
reduce the number of iterations and the processing time. 

We have used an addition uniform distribution to 
account for noise and outliers. The weight, w, of this 
distribution provides a flexible control over the method 
robustness and allows accurate CPD performance, espe- 
cially in presence of severe outliers and missing points. 

We have tested CPD on various synthetic and real 
examples and comare it to LM-ICP (in rigid case) and 
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TPS-RPM (in non-rigid case). CPD shows robust and 
accurate performance with respect to noise, outliers and 
missing points. Our method is of general interest with 
numerous computer vision applications. We provide the 
Matlab code of the CPD algorithm free for academic 
research. 
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